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ABSTRACT 

Typical flows in stellar interiors are much slower than the speed of sound. To follow the slow evolution 
of subsonic motions, various sound-proof equations are in wide use, particularly in stellar astrophysical 
fluid dynamics. These low-Mach number equations include the anelastic equations. Generally, these 
equations are valid in nearly adiabatically stratified regions like stellar convection zones, but may not 
be valid in the sub-adiabatic, stably stratified stellar radiative interiors. Understanding the coupling 
between the convection zone and the radiative interior is a problem of crucial interest and may have 
strong implications for solar and stellar dynamo theories as the interface between the two, called 
the tachocline in the Sun, plays a crucial role in many solar dynamo theories. Here we study the 
properties of gravity waves in stably-stratified atmospheres. In particular, we explore how gravity 
waves are handled in various sound-proof equations. We find that some anelastic treatments fail 
to conserve energy in stably-stratified atmospheres, instead conserving pseudo-energies that depend 
on the stratification, and we demonstrate this numerically. One anelastic equation set does conserve 
energy in all atmospheres and we provide recommendations for converting low-Mach number anelastic 
codes to this set of equations. 
Keywords: stars:interiors - Sundnterior 



1. INTRODUCTION & MOTIVATION 

In astrophysical fluid dynamics, the evolution time of 
the fluid flow is often substantially longer than the sound 
crossing time of the system. This is particularly true 
for convection deep in stellar interiors where the flows 
are very subsonic. Near the base of the solar convection 
zone the sound speed is about 220 km/s, while the con- 
vective velocities are likely of order hundreds of meters 
per second. Following the evolution of sound directly 
imposes crippling computational limits on simulations of 
such flows, as their evolution times are typically many 
convective turnover times, each of which is often several 
thousand sound times. 

So called "sound-proof" equations address this sep- 
aration of scales by beginning with the Navier- Stokes 
equations and filtering out fast, high-frequency sound 
waves while retaining compressible motions on slower 
time scales due to gravitational stratification. These 
motions include gravity waves in stably stratified re- 
gions and asymmetric convection in unstably stratified 
regions, with typically broad slow upflows and nar- 
row fast downflows. In astrophysical and geophysical 
settings, the most commonly employed "sound-proof" 
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equations are the anelastic equation s (jBatchelorl 119531 : 
lOgura fc PhTh^[l962t lGough|[l969l) . These have been 
employed in various astrophysical and geophysical 
codes to study solar convect i on an d the solar d y namo 
(e.g iGilman fe Glatzmaierl 119811: [Glatzmaierl 11984 
1985 : lClune et al.lll999| : [Miesch et aljboOOl: [Elliott et al.l 
2000: iBrun fc Toomrei 120021: iBrun et alj 120041) . stel- 
lar convection and dynamos (e.g.. iBrowning et al.1 



200 41 IBrun et al.1 



Nelson et al 



tures (e.g. 



nd d ynamos (e-g ., drowning et al.l 

20051 : IBrown"et^ 120081 I2010L 120111 : 



201 ID. the buoyan t rise of magnetic struc- 



Lantz fe Fan! I1999T) . terrestrial convection 
and the geodynamo ( e. g.. iBraginskv fe Robertsl H995I: 
Glatzmaier fe Robertsl 119961: iRoberts fe Glatzmaierl 
20001 : lOlson fe Christensenl l2006t Uones et al.l 120091: 



Jones fe Kuzanvanll2009T )" and the coupling of an unsta- 



bly stratified c onvection zone to a stably stratified region 
beneat h (e.g., 
2005bllal. [20061: 



Rogers et al.l 120031: IRogers fc Glatzmaierl 
Browning et al l 120061: IRogers et al 1120061 



20081 IRogers fc MacGregorl 120111: ferun et all I20L1 . 

Recently a significant benchmarking effort has been 
undertaken to compare the various imple mentations of 
the anelastic equations (jJones 

Formally the anelastic approximation is only valid for 
an adiabatic or nearly adiabatic atmosphere. The solar 
convection zone is nearly adiabatic but it is underlain 
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by a stably stratified radiative zone; unsurprisingly the 
anelastic equations are often extended into this region 
where their validity may break down, to study the cou- 
pling of p enetrative convection with a stably-stratified re- 
gion Ce.g.lRogers fe Glatzniaierll2006l IRogers et al.ll2008t 
iRogers fe MacGregorl2010Ll2011l : lBrun et al.ll2011h . This 
is particularly important in simulations of the solar dy- 
namo, as the stably stratified internal boundary layer 
known as the tachocline at the base of the convection 
zone is thought to play a major role in the global-scale 
dynamo. 

Fundamentally, the anelastic equations filter sound 
waves by modifying the continuity equation of the fully 
compressible Navier-Stokes equations. Questions about 
the energy conserving properties of the anelastic approx- 
imation have remained a thorny issue in the fluid dy- 
namics community, with an especially vi gorous debate 
occurring in th e atmospheric sciences (e.g., (Durran1ll98a 
iBannonl 119961 ). where these equations were originally 
derived. Likewise, there are several competing anelas- 
tic app roach es, including "co-density" for mulations (e.g., 
ILantzHl992t IBraginskv fe Roberts! Il995l hereafter the 
LBR equations) and their different properties are un- 
clear. 

An alternate approach to sound-proofing the Navier- 
Stokes equations are the pseudo-incompressible equa- 
tions, where the pressure rather than continuity e quation 
is mo dified. These equations were proposed in (jDurranl 
l!989T l and have recently been adopted in the astro- 
physic al fluid dynamics comm unity (e.g.. lAhngren et al.1 
I2006allbl ; iZingale et a 11 l2009D and se e particular use in 
the MAESTRO code (jNonaka et al.l 1201(1 . The prop- 
erties of gravity waves and stable-layer dynamics in the 
pseudo-incompressible equations have been explored ex- 
tensively in the atmospheric sciences community, with 
several comp arisons agains t the properties of the anelas- 
tic equations (iDurranl [19891 120081 INance fc Durranll 19941: 



lAchatz et al.ll2010HKlein et al.ll2010ft . We reserve further 
discussion of gravity waves in this set of equations for a 
later paper. 

Here we explore three implementations of the anelastic 
equations, one used in the anelastic spherical harmonic 
(ASH) code, and two different implementations of the 
"co-density" formulation (LBR equations). These equa- 
tion sets are detailed in Section [2j We show that the 
anelastic equations based directly on the Navier-Stokes 
equations (anelastic Navier-Stokes, or ANS equations) 
behave incorrectly in stably stratified region. First we 
analytically study wave motions in an isothermal atmo- 
sphere in Section EH We find that these equations do 
not conserve energy and instead conserve an entropy- 
weighted "pseudo-energy" (Section |4|). We find however 
that the LBR equations do behave correctly for strongly 
stratified regions, conserving energy and reproducing the 
results obtained from the full compressible Euler equa- 
tions. This is surprising, as the LBR equations make 
further assumptions of adiabaticity beyond those con- 
tained in the basic ANS equations, but these assump- 
tions lie at the heart of the energy-conserving proper- 
ties. As a consequence, adjustments to the LBR equa- 
tions to more correctly capture the sub-adiabatic strat- 
ification can have profound consequences, introducing 
a co mpletely different form of en ergy non-conservation 
(e.g., IRogers fc Gla tzmaier 2005b! an d hereafter the RG 



equations). We explore the behavior of these differing 
equations further in bounded atmospheres and spherical 
geometries in Section [5] and perform numerical simula- 
tions that show the difference between the normal ANS 
equations and the LBR equations. The implications of 
these findings for simulations of solar convection is dis- 
cussed in Section [51 which also give suggestions for im- 
proving anelastic treatments of stably-stratified regions. 
The reader who is primarily interested in implementing 
energy-conserving anelastic equations should read Sec- 
tions El M and El 

2. MODEL EQUATIONS 

2.1. Fully compressible Euler equations 

For the purposes of this paper, the most general equa- 
tions for fluid dynamics in the solar interior are the fully 
compressible Navier-Stokes equations. When viscosity is 
neglected, as we do here, these are known as the fully 
compressible Euler equations (FC equations). The equa- 
tions of continuity and momentum are 



dp 

— + u- Vp=-pV u, 

at 

du \ 

— + u-Vu)=-VP + P g, 



(1) 
(2) 



with gravitational acceleration g = —gr. For an ideal 

gas, 

P = IZpT = (7 - l)p£, (3) 

with £ the specific internal energy and 7 = cpjcy = 5/3 
is the ratio of specific heats. Here cp = 7/(7 — 
is the specific heat at constant pressure. The evolution 
equations for temperature and pressure are 

^-+u-VT = -( 7 -l)TV-u, (4) 
at 

dP 

— +u ■ VP=-7PV • u, (5) 
at 

where thermal conduction and other sources and sinks of 
energy are neglected. 

Although equations ([T]-[5]) form a complete system, it 
will be useful during our discussion of the anelastic equa- 
tions to rewrite these in terms of entropy S. Equa- 
tions @ and ([5j) can be combined with an equation of 
state linking the thermodynamic properties 

— = -d\nP-d\np= -dlnT- — -d\np (6) 
cp 7 7 7 

into an equation for the evolution of entropy fluctuations 

m +u - vs = °- W 

We now specialize to the case of a hydrostatically bal- 
anced, stratified atmosphere with background density 
stratification po, pressure Pq, temperature Tq and en- 
tropy Sq that only vary with radius, with 



VP = Pog- 



(8) 



We define fluctuating quantities, denoted with subscript 
1, by subtracting the time- independent hydrostatic at- 
mosphere making no assumptions about relative ampli- 
tudes, with e.g., Pi = P — Po(r), thus these equations 
are fully nonlinear. 
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2.2. Anelastic models and fully compressible Euler 
equations in standard form 

All anelastic approximations employ a continuity equa- 
tion of the form 

V • ( Po u) = 0. (9) 

Equation @ derives from the assumption that the den- 
sity fluctuations are small 



pi = p - po < po- 



rn 



(11) 



In this case, the fluctuating density is given by the lin- 
earized equation of state, 

Pi = 1 P_ _ Si = Pi _Tl 

Po J Po c P P T ' 
and though using a linear equation of state is not strictly 
required, we find it a clarifying simplification for the cur- 
rent discussion. 

We consider equation (O to be the defining character- 
istic of anelastic models. There exist however a variety of 
different treatments for the momentum and energy equa- 
tions in the anelastic literature. In the following subsec- 
tions we will consider three common formulations. The 
different notation and different thermodynamics used in 
the various anelastic treatments leads to some confusion. 
To remedy this, we reproduce each set of models under 
as consistent a notation as possible. Practical numerical 
or computational differences can arise when solving dif- 
ferent transformations of the same fundamental model, 
but these issues lie beyond our current scope. Therefore, 
we consider two models identical if one can bring them 
into the same form by legitimate mathematical transfor- 
mation, i.e., without approximation. 

For comparison with the anelastic equations we first 
write the FC equations in standard form. With a lin- 
earized equation of state (jllj) . we rewrite the buoyancy 
term involving pressure fluctuations in the following fash- 
ion 

Pi P a VP Pi 



-g = 

7P0 Po 7"o 



Po 



{ So 

cp 



+ Vlnpo 



(12) 



where we have used equations ([SJ and ([5]). We now in- 
troduce the reduced or kinematic pressure w with 

w=— . (13) 

The fully compressible Euler equations, with an entropy 
based energy equation and reduced pressure w, are 

Vp = -poV-M, (14) 

at 

du „ „{S \ Si , x 
— +u-Vu = -Vvj + tuV — -g, 15) 

at \c P c P 



dSi 
dt 



u- VSi = -u ■ VS . 



(16) 



These equations linearize the thermodynamic variables 
(eq. Hip but are nonlinear in the velocities and are the 
counterparts of the anelastic equations that we now turn 
to; we do not solve these equations (TI4T[T6"]) but include 
them for illustrative purposes. 

2.3. ANS Anelastic equations 

In many anelastic equations the momentum 
equati on is the same as in t he FC equations 
(e.g., iGilman fc Glatzmaierl 119811: iDrew et all 119951: 



iClune et alj|1999t iBrun et all 120041 ). We thus refer to 
these equations as the anelastic Navier-Stokes (ANS) 
equations. In the ANS equations, the momentum 
equation is 



I ^ + («• V)it 



- VPi + pig, 



(17) 



which with equation (jllj) can be transformed into the 
same form as equation (|15p . with 

^.+u-Vu = -Vw + wv(^)-^g. (18) 
dt \cpj c P 

The ANS momentum equation (|18p can be written in an 
alternative form by combining w terms to yield 

^ + u ■ Vu = -e«Vcp) V (W^o/cpA - ^g, (19) 
dt V / cp 

which will be useful for our analysis in Section [4] As a 
notational issue, the stratification term interacting with 
w in equation (|19[) takes the same form as a potential 
temperature, as is traditionally used in studies of geo- 
physical flows in the atmosphere and ocean with 



e = e^ c ^ 



p 



1/7 



Po 



(20) 



(21) 



and, with the linearized equation of state (|11[) . 
Oi = Si 

60 c P 

In terms of 0, the ANS momentum equation is 

^L + u-vu = -e v(vue Q 1 )-^g. (22) 

Ot fo 

Neglecting diffusion and sources of energy, the energy 
equation is the same as in the FC equations (eq. [7} , with 
a background entropy gradient 

-^+u-VS 1 = -u-VS . (23) 

Combined with the anelastic continuity equation @, 
equations (fT8|) and (|23|) constitute a full set of equations 
for anelastic motions. 



2.4. LBR Anelastic equations 

In the "co-density" equations or 
Braginsky-Roberts eq uation s fe.g.. lLantzl 
Braginskv fc Roberts! 119951 : lLantz fc Fan! 
Jones et al. I I2009L and hereafter LBR 



Lantz- 



1992 



1999 



equations), 
momentum 



the n7V(5'o/cp) term is dropped and the 

equation becomes 

du _ _ iSi . „. 

— +u-Vu= -Vm g, 24 

at c P 

or, in terms of potential temperatures, 

— +«■ Vu= Vw-— g. 25 

at Bo 
As in the ANS equations, the energy equation (|2"3"j) and 
the continuity equation Q complete the full set of equa- 
tions. 

The LBR momentum equation (|24p is derived from 
the full Euler momentum equation ^ by assuming 
VSo w 0, as for a nearly adiabatic state. Despite 
this assumption, we find that the LBR equations per- 
form well when V5o 7^ while the ANS equations per - 
for m poorly in that same limit Though lLantzl (|1992ft 
and iBraginskv fc: Roberts! (|1995l ) are typically credited 
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with independently deriving the LBR equations, these 
equations bear striki ng similarities to the Lipps-Hemler 
anelastic equations (jLipps fc Hemlerl 119821 119851 iLippsI 
ll990T l , who were possibly the first to recognize the impor- 
tance of introducing a reduced pressure and neglecting 
the interactions between fluctuating pressure and strat- 
ification. They likewise recognized that gravity waves 
derived from their anelastic equations conserved energy. 

2.5. RG anelastic equations 

iRogers fc Glatzmaierl (|2005bf ) use a different set of 
anelastic equations (hereafter the RG equations). As 
above, neglecting viscosity and sources of heat, their 
equations are the momentum equation and a tempera- 
ture based energy equation 

T\ 

tuV In Tq g, 



— — h u ■ Vu = — Vt 



du 
~dt 

m_ 

~~dt 



u ■ VTi 



To' 
l)(To 



(26) 



it - VT - (7- l)('ib + iiJV • u 

(27) 

With the anelastic continuity equation (|9]), these consti- 
tute a full set of equations for anelastic motions. 
Equation (|26|) can be equivalently written 

^+u- Vu = -T V(w/T )-^g, (28) 
at j 
a form that will be useful in Section [4] By combining 
the equation of state ([6]) with the anelastic continuity 
equation we can cast the energy equation in terms 
of entropy with 



d_ 

dt 



T 



7 1 + 



With a linearized equation of state, this takes the same 
form as the entropy equation (|16p . but with an extra 
factor of 7 multiplying the background entropy gradient. 

The right hand sides of the momentum equations for 
these four systems of equations are summarized in Ta- 
ble ni 

3. LINEAR WAVES IN AN ISOTHERMAL 
ATMOSPHERE 

A plane-parallel isothermal atmosphere gives an an- 
alytically tractable background for computing cigenfrc- 
quencies and modes for linear gravity and/or acoustic 
waves. Computing these simple solutions helps elucidate 
the differences between the various anelastic treatments. 
Defining the velocity in terms of the vector displacement, 

H (30) 



u 



dt 



allows simple integration of the linear thermodynamic 
equations, 

Pi/po = -€-Vlnpo-V.€, (31) 

Pi/lPo = -£- VlnP 1/7 - V-£, (32) 

Si = -e-VS . (33) 

For wavelike perturbations in an atmosphere of infinite 
extent, we can assume without loss of generality that 

(£, pi,P\,Si) cx f(Kr) exp (iiut — imx) (34) 

where x is the horizontal coordinate, m is the horizontal 
wave number, and the vertical dependence on r has been 
left in general form with wavenumber K. 



Table 1 

Systems of equations 



System 


RHS momentum equation 


cq 


FC 


-Vw + wV(S /cj>) - (S 1 /c P )g 




ANS 


-Vro + roV (So/cp) - (Si/c P )g 




LBR 


-Vc7- (Si/c P )g 




RG 


- Vro + ro V In To - (Ti /T )g 




RHS wave momentum 


FC 


-Vvj + roV (So/cp) + (£ ■ V)(S /c P )g 




ANS 


-Vm + roV (S /cp) + (£ ■ V)(5o/c P )g 




LBR 


-Vti7+(€- V)(S /cp)s 




RG 


- Vro + roVlnTo + T (g ■ V)(5 /c P )g 




RHS momentum (for Section [4} 


ANS 


_ e S /cp v (roe- s o/ c p) - (Si/cp) g 




LBR 


-Vict - (Si/c P )g 




RG 


-T V( ro /To)-(T 1 /To)g 


(J26jl 



Note. — In all systems of equations, to = Pi/ pa- 
The fully compressible equations use continuity equation GJ 
while anelastic systems use equation $9§. In the wave mo- 
mentum equations, £ is the displacement vector as defined 
in eo 1301 

In a hydrostatically balanced isothermal atmosphere 

9 l_ _ 75 

4 



V r In Pq = V r In pa = — 



TZT ( 



ii 



H 



(35) 



where H is the pressure or density scale height, V r is the 
vertical derivative, and 



c 2 ^ 



(36) 



is the (constant) sound speed. The Brunt- Vaisala fre- 
quency is 



N 



-gV 



(7 ~ 1) 9 
7 H 



(37) 



3.1. Fully compressible equations 

The solution for the fully compressible equations is 
well known an d can be found in several textbooks (e.g., 
lLighthilj||1978[ ) . We begin with the linearized momentum 
equation for waves 



- Po u 2 £, = -VP l +p 1 g, 



(38) 



and solve for eigenfrcquencies using equations (I30l) - (|32p 
and (|54")) . Taking the vertical eigenfunction f(Kr) = 
exp (— iKr), the dispersion relationship for waves in an 
isothermally-stratificd atmosphere is 



\K 



iKH~ 



m 2 N 2 . 



(39) 



It is well known that the fully compressible Euler equa- 
tions conserve energy. Their frequencies w must thus be 
purely real with no imaginary component (see Section|4]), 
yet equation (|39j) has an imaginary component. Taking 
a complex vertical wavenumber 

K = k + l ^H (40) 
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with real component k resolves this. The vertical eigen- 
function follows 



f(Kr) 



exp (—iKr) = exp ( — — ) exp {—ikr). 
V 2H I 



(41) 



All waves in this atmosphere share the properties that 
their eigenfunctions grow with height, their momentum 
density decreases with height 

p u oc exp (~^), (42) 

while their kinetic energy p^u 2 is constant with height. 
These eigenfunctions are orthogonal with respect to the 
density weight 



p f(Kr)f(K'r)*dr = S(k-k'), 

with 5 here the Dirac delta. 
The final dispersion relationship with w 2 real is 



1 



AH 2 



= m 2 N 2 . 



(43) 



(44) 



with N 2 given by equation (|37j) . The quadratic nature of 
equation (|44[) in w 2 provides for two distinct branches of 
acoustic and gravity waves. In the high frequency limit 
oj 2 » A 2 , 



2 



1 



AH 2 



(45) 



representing the propagation of pure sound waves in an 
atmosphere with an acoustic cutoff frequency c| / AH 2 . In 
the low frequency limit, we obtain pure internal gravity 
waves with 

,2 



2 



The full solution for 



k 2 

j 2 follows 



l 

TW- 



-N 



J sw 



1 ± 



l 1 _^GW 



(46) 



(47) 



J sw 



with the positive root corresponding to the sound waves 
while the negative root corresponds to the internal grav- 
ity waves. 

3.2. A NS gravity waves 

We begin our analysis of the anelastic systems with the 
ANS equations. For linear waves, the continuity, momen- 
tum and energy equations are 



Vlnpo, 
-Vtn + wV 



So \ Si 

— g, 

CP J C P 



Si = -£ • VS 



(48) 
(49) 
(50) 



Combining the vertical momentum equation (|49|) and the 
energy equation (|50|) for linearized waves, we obtain 



UJ 2 ^ r 



V r (vu) + vuV r (So/cp) 



(51) 



where we have also used equation (|37|) . We obtain zu 
by taking the horizontal divergence of the momentum 
equation 

lu 2 V_l • £ = V\vj. (52) 



The dispersion relationship for linear waves is 

2 



K 



.2 7 - 1 
'' 2-/H 



1 



4 7 2 iJ 2 



m 2 A 2 . (53) 



Once again, lo 2 has an imaginary component. As with 
the fully compressible equations, we can try to absorb 
this imaginary component within the vertical eigenfunc- 
tion, which leads to a vertical wave number 

1 2 7 - 1 

(54) 



K = k + i 
and vertical eigenfunction 



2H 



1 



f{Kr) = exp (—iKr) = exp 



2 7 - 1 



2H 



with time dependence 

i 2 + k 2 



exp (—ikr), 
(55) 



2 

W ANS 



1 



= m 2 N 2 . 



^ 4 7 2 ;PJ-'v • ^ 

A serious problem lurks within these choices however, as 
the momentum and kinetic energy densities scale as 



~2^H 



PqU occxp 



1 r 

7 — 1 r 
7 H 



(57) 
(58) 



For adiabatic motions in an ideal gas, 7 = 5/3, and 
the kinetic energy of the waves grows exponentially with 
height. 

Alternatively, we can use the eigenfunctions from the 

fully compressible equations, in equations (|41ti40|) . which 

leads to the correct far-field behavior for momentum and 

energy, but leads to a dispersion relationship of 

1 2 — t 1 

• ikg(-i-l) =m 2 N 2 . 



W ANS-I 



AH 2 



7 



(59) 

There is now an imaginary component to lo 2 and anelas- 
tic gravity waves in an infinite isothermal atmosphere 
can have spurious growing (or decaying) modes. As we 
will see in Section[4] this bizarre behavior reflects the fact 
that the ANS equations do not conserve energy. Further, 
as we will see in Sections 01 13 the fact that these spurious 
modes have not been detected in simulations previously 
is likely related to the presence of a conserved pseudo- 
energy (i.e., a differently weighted quadratic integral of 
the fluctuating velocities and entropies). 

3.3. LBR gravity waves 

Finding linear eigenfrequencies in the LBR equations 
amounts to the same procedure as in Section 13.21 Now 
however the tjjV (Sq / cp) term is missing from the verti- 
cal momentum equation, and equation (|51[) becomes 

w 2 £ r - V r n7 = iV 2 £ r . (60) 

This readily yields the following dispersion relationship 

1 \ 2 1 



K - i- 



2H 



AH 2 



m 2 N 2 



(61) 



Requiring that K = k+i/2H is clearly the natural choice 
for obtaining real frequencies. By employing the vertical 
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Figure 1. Dispersion relationships for gravity waves in an isother- 
mal atmosphere of infinite extent with 1/H = 10, for the fun- 
damental mode k = 1 and with horizontal wavenumber m. (a) 
Frequencies for each set of equations in the low frequency limit 
(ui 2 /N 2 -C 1). Shown are the gravity wave branch of the exact so- 
lutions for the fully com pressible Eulcr equations (black, labelled 
FC, and given by eg I47|l . Also shown arc the dispersion relation 
for the ANS equations (blue, dot-dashed), the LBR equations (red, 
solid), and the RG equations (green, dashed) with each line labeled. 
In this regime the LBR equations and exact solutions to the Eu- 
ler (FC) equations are in good agreement, while the ANS and RG 
equations obtain frequencies that are too high. The corresponding 
dispersion relationships arc given in Table[2] (6) Full frequency do- 
main. At large m the ANS, LBR and FC equations converge to the 
Brunt- Vaisala frequency N, while the RG equations are too large 
by a factor of ^7. Here we also show t he s ound wave branch (black, 
labelled "acoustic" , and given by eq 1470 of the exact solution to 
the full Euler equations. 



cigenfunctions in equations (|41M0j) we get 

1 



2 

W LBR 



mf + k + 



4H 2 



m 2 N 2 , 



(62) 



which is the same as equation (|itj| . Adiabatically prop- 
agating gravity waves solved with the LBR equations in 
an infinite isothermal atmosphere behave like the low 
frequency branch of the fully compressible equations in 
both their time dependence and their vertical structure. 

3.4. RG gravity waves 

Next we look at the propagation of gravity waves 
within the RG equations. In an isothermal atmosphere 
the coupling between w and the background stratifi- 
cation disappears. With the anclastic continuity equa- 



, the linear RG wave equations are 

2f xy Ti 

-lu £ = -Vm- —g, 
Jo 



(63) 
(64) 



-i = - 7 £-V(So/c P ). 
. 

Combining the vertical momentum and energy equations 
yields 

0J 2 & - V r nj = jN 2 £ r (65) 



This leads to a dispersion relationship of 

^ 2 

K - i 



1 
2H 



1 



AH 2 



(66) 



As previously, the vertical eigenfunctions in equa- 
tions (|4HI40[) are the clear choice and lead to a final dis- 
persion relationship of 



J RG 



m 2 + k 2 + 



1 



AH 2 



m 2 <yN 2 , 



(67) 



which is the same as equation (|46| except for the factor 
of 7 multiplying N 2 . 

While the functional form of the frequencies given in 
equation (|67)l are correct up to a factor of ^7 w 1.29 
for 7 = 5/3, and while the vertical structure of the 
eigenfunction matches with the fully compressible case, 
we note that this is a special case brought about by 
V In To = in an isothermal atmosphere. In more gen- 
eral atmospheres, an extra term would exist in equa- 
tion (f63|) of the form w V In To and we would be faced by 
the same problems with energy conservation and growth 
that we found in Section [3~2l for the ANS equations. We 
will see this in Section 1431 

We summarize the properties of gravity waves for all 
four systems of equations in an isothermal atmosphere of 
infinite extent in Table [2] and plot them for waves with 
kH = 1/10 in Figured] In the low-frequency regime 
(Figure [TJi), the gravity wave branch of the exact solu- 
tions to the Euler equations (labelled FC, and given by 
cq RTF]) matches the dispersion relationship of the LBR 
equations closely, while the frequencies of gravity waves 
in the ANS and RG equations are too large. As the 
horizontal wavenumber m increases, the LBR dispersion 
relationship begins to diverge from the exact results. At 
still larger wavenumber m, both the ANS and LBR dis- 
persion relationships return to agreement with the exact 
Euler solutions (Figure QJ). At all wavenumbers, the 
frequencies from the RG equations are a factor of ^7 
larger than those obtained from the LBR equations and 
thus exceed the Brunt- Vaisala frequency N at large m. 
Higher order radial modes show similar behavior, though 
the relative differences between the LBR and FC disper- 
sion relationships decreases as k increases. 

4. CONSERVATION OF ENERGY AND 
PSEUDO-ENERGY 

An important theme of this paper revolves around en- 
ergy budgets in different approximations to the full Euler 
equations. The curious discrepancies found in isothermal 
atmospheres in Section [3] hint at deeper issues in these 
approximated equation sets. In this section, we find that 
those issues are associated with energy conservation and 
its violation. Here we consider general atmospheres, with 
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Table 2 

Infinite isothermal atmosphere 



System 



eq 



FC 
ANS 
ANS-I 
LBR 
RG 



m 2 N 2 
■m 2 N 2 



l 2 — v 



^2=2+^(7-1)1 m 2 N 2 



m 

2 . 



1 

AH' 



m N 



m 2 ^N 2 



l62ll 



Note. — Horizontal wavenumbers m and vertical wave num- 
bers k are real quantities. In all systems of equations except ANS, 
we ha ve take n the radial eigenfunctions corresponding to equa- 
tions l|4H440jl which remain finite in the far-fiel d limit. In the 
ANS equations, we instead take eigenfunctions I54H55H . which 
leads to real lo but divergent behavior in the far-field limit. For 
the FC equations we here show only ^q w (the low-frequency 
limit ); t he full Euler dispersion relationship is given in equa- 
tion B7t . 



Table 3 

Energies and pseudo-energies 



System 



weight 



eq IA RZ CZ 



FC po 
ANS exp(-5o/c p )po 
LBR po 



i 75i 
iggl 
96l 



RG 



'o Po 



Y Y Y 
N N Y 

Y Y Y 

Y N N 



Note. — Weight required for sclf-adjointncss and 
hence energy or pseudo-energy conservation in each 
system of equations, with reference to where the con- 
servation properties are shown in the text. Systems 
with weights other than po will not always conserve 
energy. Included arc qualitative estimates of whether 
each set of equations is likely to conserve energy in an 
isothermal atmosphere (IA), in the stably stratified 
solar radiative zone (RZ), and in the nearly adiabat- 
ically stratified solar convection zone (CZ). 

the isothermal atmospheres of Scction[3]bcing a subset of 
these results. For each set of equations, beginning with 
the full Euler equations and proceeding with each anelas- 
tic equation set in turn, we derive the energy conserva- 
tion properties for arbitrary nonlinear motions. We then 
consider the energy conserving properties of linearized 
motions including wave-like perturbations. We find that 
some equation sets (FC and LBR) conserve energy and 
behave as expected. We find that the other anelastic 
equation sets (ANS and RG) do not conserve energy and 
instead conserve a stratification- weighted pseudo-energy, 
which leads to some surprising and paradoxical results 
for wave-like motions. The key results of this section are 
summarized in Table [3] 

4.1. Euler Energy Balance 

We begin by considering the Euler equations. The 
main results of this subsecti on are well known in the lit- 
erature (e.g.. iLighthilll [19781 ). Namely, in the fully com- 
pressible Euler equations, energy is conserved by wavclike 
motions and the temporal frequencies oj are purely real. 
However, for the purposes of comparison with anelastic 
models, we note that the fully nonlinear equations (Q]) 



((71) contain a statement of conservation of energy. Con- 
tracting equation ([2]) with velocity u and assuming that 
gravity is given by a potential function gives 
r) E 

^ + V-[u(E + P)} = 0, (68) 



where 



dt 



E = 



P\u\ 2 
2 

-V$, 



P 



7-1 



(69) 

g = -V$, (70) 

with 4> the gravitational potential. The fully compress- 
ible Euler equations conserve energy for arbitrary (non- 
linear) motions. 

For the linearized version of the Euler equations, we 
may go a step further. For a system in hydrostatic bal- 
ance (eq. [5]), we write equation Q38p in terms perturbed 
pressure P\ and entropy Si as 



POgi = — VPl 



Pi q 2 On Si VSn , , 

^VPo + ^^-^. (71) 



c P c P 



We introduce an arbitrary vector that is related to the 
displacement vector £ (eq 150)) and guided by equations 
((321) and ((33]) define P[ and S[ as 



P[ = -t'- VP - 7 PoV-£', 



(72) 
(73) 



Contracting equation (|7ip with arbitrary and using 
equations ([72"]) and ([73]) gives 

at z jPo N z cp cp 

(74) 

Wc may derive a number of different results from equa- 
tion (j74[) . First we consider velocity perturbations and 
take = dti = u (thus S[ = d t Si and P{ = d t Pi). This 
choice gives the local conservation of energy for linear 
perturbations 



2 / o \ 2 
9 Po I Si 




Pi 2 
27P0 



V • (uPi) = 0. 

(75) 

Integrating equation ([75]) over a volume V with (u ■ 
n) Pi = on the boundary dV, gives 

^(K + U) = 0, (76) 

where the kinetic and potential energies are given respec- 
tively by 

(77) 



K=^j^p Q \u\ 2 d 3 x, 



U 



1 f fg 2 po fsA 2 , p 2 



2 Jy \ N 2 V On 



+ I d 3 *. (78) 
7P0 / 



Linear perturbations also conserve energy in the fully 
compressible Euler equations. 

Choosing instead that = £ (with this choice, S[ = Si 
and P'i = Pi) in equation (]74[) . integrating over volume 
V with (u ■ fi) Pi =0 on boundary dV , and averaging 
over time gives a version of energy equipartition for lin- 
ear perturbations, where the time average of the kinetic 
energy equals the time average of the potential energy. 
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Rather than considering energy conservation, we now 
consider the timc-dcpcndcncc of linearized displace- 
ments 



€ = r 



(79) 



where £* represents the complex conjugate of displace- 
ment £. Here, S[ gives the characteristic entropy pertur- 
bation Si associated with a displacement of amplitude £ 
and likewise with P[ and pressure perturbation P\. In- 
tegrating over the same volume, V, gives 



Jv Jv 



9 2 Po \Si\' 
N 2 cl 



d'x = 0. 



(80) 

All of the integrals in equation ([80]) are strictly real and 
positive definite, which implies that the squared tempo- 
ral frequencies must also be real 

3 (w 2 ) = 0. (81) 



Equation (|81[) states that while instability may or may 
not exist, the system must transition from purely oscil- 
lating (TV 2 > 0) to purely growing behavior (N < 0). 
Neither growing nor damped waves exist in the fully com- 
pressible Eulcr equations. 

Lastly, one may show that displacements with differ- 
ent frequencies are orthogonal with respect to the energy 
inner product, 

(82) 



where 5 u ' lU1 is here the Kroncckcr delta. Together, equa- 
tion l]82p and equation (|74|) imply that the right-hand 
side of the linear perturbation equation (|71|) is self- 
adjoint with respect to this energy inner product. There- 
fore, the condition equation (j8"Tj) unlimitedly stems from 
both a particular dynamical equation, and an appropriate 
inner product. If equation ()82[) is altered, which amounts 
to a different spatial weighting of the solutions, then 
equation (|8ip may not hold, and the time dependence 
of the solution may acquire spurious growth or decay. 

The above four results that derive from integrating 
over the volume V hinge on the condition that 



(it • n) Pi = 



(83) 



on the boundary of V with n the unit normal vector. 
This condition is not a mere technical triviality, as equa- 
tion (|83|) causes the divergence term in equation (|74| to 
vanish. If V is a bounded domain, or is periodic in the 
horizontal direction and bounded in the vertical direc- 
tion, then we may easily satisfy equation (|83|) by requir- 
ing (u • n) — by itself (e.g., impenetrable boundaries). 
For the travelling waves we considered in Section |3] the 
product (u-n) P is itself periodic and integrates to zero, 
since \u\ ~ p , and |Pi| ~ p for large and small at- 
mospheric heights. As we will see in the following subsec- 
tions, the far-field behavior of travelling waves controls 
the stability properties of different anelastic models. 

4.2. ANS Energy Balance 

For comparison with the total energy equation (|68[) for 
the Euler system, we now derive an equivalent energy 
balance for the anelastic models, beginning with the ANS 



equations. Contracting equation (|18[) with u and using 
the anelastic continuity equation (|9]) gives in basic form 

_ . ,__ Si _ f So 

— + V-[u{K + p w)\+p u-g — = zjp u- V — 

(84) 

with kinetic ene rgy density K = pa\u\ 2 /2. Using the 
relationship IA2I in Appendix [K\ we put the left hand 
side of equation (|84[) into conservative form 

dE f S \ 

— + V-[u{E + p m) ] = vopou ■ V i^j , (85) 

where E and w are given by equations (IA4[) and (IA5|) 
respectively. We cannot however transform the right- 
hand side into conservative form unless 

lim — / / wpou • VS d 3 xdt = 0. (86) 

This condition is not true in general and, simply stated, 
arbitrary (nonlinear) motions in the ANS equations do 
not conserve energy. Condition (|56"|) is satisfied for 
adiabatically-stratified atmospheres, where VSo = 0, 
and in those systems the ANS equations do conserve en- 
ergy. 

We turn now to linearized motions to learn more about 
the strange behavior found in Section [3. 2 1 by considering 
the equivalent of equation (|74p for the ANS model equa- 
tions. Contracting the linear momentum equation with 
an arbitrary but here satisfying V • (poO = (again, 
£' could be either £ or <9 t £ = u), produces 



dt 2 N 2 cl 



-V-{ P oi'w)^mpoi'-V 



So 



(87) 

If we integrate equation (|87|) over the entire volume, V , 
then the right-hand side refuses to vanish: even linearized 
motions do not conserve energy in the ANS equations. 

The non-vanishing right-hand side of equation (|8T[) 
would also appear to imply that the squared frequencies 
ui 2 are not strictly real. On the surface, the asymmet- 
ric nature of equation (|87|) would appear to imply non- 
sclf-adjointncss of the linear equations and hence spuri- 
ously growing modes. This is consistent with what we 
found for our analysis in an infinite isothermal atmo- 
sphere (Sec. 13. 2p ; as we found there, a correction to the 
spatial structure counteracts this effect and regains real 
eigenvalues for the linear equations at the cost of modes 
which grow in spatial height. In the literature of anelastic 
simulations however, no mention appears of these spuri- 
ously gro wing gravity waves, and a paradox seems appar - 
ent (e.g- lRogers fc Glatzmaierll2005bl : iBrun et alJl20111 ). 

The paradox of spurious growth is remedied by the 
following transformation of equation (|87|l . 

, ( t , &t , 9 2 S[Si 

where we define the scaled pseudo-density 

Po^Poe- So/c ", 

which reduces to the actual background density in the 
case of adiabatic stratification. 

Though energy is not conserved for nonlinear dynam- 
ics, nor for linear waves, equation (|88|) implies that the 



V-OSo£V) = 0, 



(89) 
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following pscudo-cncrgy is conserved 

'Si 



2 



Po 



u 



9 

N 2 



d 3 x, 



(90) 



i.e., dtE = 0, for linear perturbations. If the pertur- 
bations are nonlinear then the rescaling of the density 
fails since the advection of kinetic energy is not an exact 
divergence in terms of this pseudo-density. 

As in the compressible case (eq. I50"|) , one may use equa- 
tions (ED) & JEIl) to show that 



/ p m 2 d 3 x 

Jv 



f n ^ 


Si 


l v Po W 2 


Cp 



d*x = 



f N 2 

/ po— II -g\ 2 d 3 x, 

Jv 9 



(91) 



which implies that 3 (w 2 ) = even if energy is not con- 
served. This indicates that the conservation of a pseudo- 
energy resolves the paradox of spurious growth and leads 
to purely real squared temporal frequencies uj 2 . We be- 
lieve that this explains why this phenomena of pseudo- 
energy conservation and energy violation has been pre- 
viously missed in the literature. 

Equations ([88]) & ([91]) imply that the linearized ANS 
equations are self-adjoint under the pseudo-energy in- 
ner product, and that eigenfunctions with different fre- 
quency arc orthogonal with respect to this pseudo- 
density weighted norm 



PoH'-td 3 x 



(92) 



The difference between equations (|82j) & (J92J) imply that 
external forcings and initial conditions project onto dif- 
ferent frequencies and basis vectors differently in the 
ANS equations than in the FC equations. In particular, 
the eigenfunctions of pscudo-energy-conserving waves in 
the ANS equations are different than the eigenfunctions 
given by energy-conserving motions (e.g., the FC equa- 
tions). In strongly stably-stratified atmospheres, these 
differences may be dramatic, as we will encounter in Sec- 
tion [5] 

4.3. LBR Energy Balance 

Unlike the ANS equations, the LBR equations show no 
problems with energy conservation. Contracting equa- 
tion (|24l) with u and using the anelastic continuity equa- 
tion §§§ gives in basic form 

■■ V.[u(JT + pdw)]+Po«-fl— =0. (93) 

c P 



at 



Using the relationship IA2I in Appendix [A] we put equa- 
tion (|93p into conservative form 
f)E 

— + V-[u(E + p Q w)} = 0, (94) 

where E and w are given by equations (|A4[) and (|A5j) 
respectively. If we integrate this over a bounded volume 
V (where as in eq [55J u ■ n = 0) then the divergence 
terms vanish and arbitrary (nonlinear) motions in the 
LBR equations obey an energy conservation law. 

For linear perturbations and for nonlinear perturba- 
tions in certain atmospheres (including adiabatic and 



isothermal atmospheres) , the LBR equations conserve an 
alternative total energy 



*-5 



Po 



u \ 2 + ±^Lst)d s x, 

c p aoo 



(95) 



i.e. d t E = 0, as detailed in Appendix [A] 
For linear perturbations we may furthermore write 



Po 



* " dt 2 



g 2 S[Si 
~N 2 ~cl~ 



V-(p £V>=0. (96) 



This implies self-adjointness of system under the energy 
inner product, and also that 



d 3 a; = 
N 2 



( 9 2 


Si 




Cp 



Po-T 

v 9 



t-g\'d 6 x, (97) 



whence it follows that 3 (w 2 ) = 0. Linear motions con- 
serve energy in the LBR equations and wavelike motions 
have real squared temporal frequencies to 2 . 

4.4. RG Energy Balance 

For the RG equations, using similar transformations 
as in Sections 14.1114.31 we obtain the following nonlinear 
energy balance 



OE 

— + V-[u(E + Po w)] 



where 



E 



Po 



2 



7$ 



ropou • VlnT , (98) 
Ti 



To 



w = w-— I $(r)dS (r) 

Cp J a 



(99) 
(100) 



It is not possible to cast the right-hand side of equa- 
tion ([55} into conservative form except in the specialized 
case of isothermal atmospheres where VlnXb = 0. This 
contrasts with the ANS equations, which can only be 
written in conservative form in adiabatic atmospheres. 
Thus the RG equations do not conserve energy for either 
anelastic convection or gravity wave dynamics in arbi- 
trary atmospheres. 

Linear perturbations to these equations do nevertheless 
obey a pseudo-density weighted sclf-adjointncss 



d 2 £ g 2 T[Ti 



V-(/3 £V) = 0, 



where 



— = - 7 £ . V — 

Jo V c v 



and the pseudo-density becomes 

Po 

Po = 



(101) 
(102) 

(103) 
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Furthermore, as in Sections I4.1H4.31 we find that 



^ / Po\Z\ 2 d 3 x 
Jv 



f * 9 2 




l v P %N 2 


T 



N 2 
'V 9 



Z-g\ 2 d d x. (104) 



Equation (|104p implies that 9 (to 2 ) =0 even if energy is 
not conserved. As in the ANS equations, waves in the 
RG equations have real squared temporal frequencies w 2 
in volumes where po (£ • n)uj = on the domain bound- 
aries, but the eigenf unctions and energies are weighted 
by pseudo-density (|103[) . 

Equation (|104[) implies that the stability boundary 
for the fully compressible and other anelastic models, 
N 2 = 0, remains unaltered in the RG equations in spite 
of energy non-conservation. Two problems do however 
still remain. The first is that an extra factor of 7 appears 
in the last integral of equation (|104[) . As we found for 
waves in an isothermal atmosphere, this leads to frequen- 
cies that are too high. The second more serious issue is 
that energy is not conserved unless the background at- 
mosphere is isothermal. In particular, both linear and 
nonlinear motions within adiabatically-stratified atmo- 
spheres will not conserve energy. 

5. BOUNDED ATMOSPHERES AND 
IMPLEMENTATION IN SPHERICAL SYSTEMS 

We now turn to considering gravity waves in a spher- 
ical shell. In this geometry divergence at infinity is no 
longer a problem. We will find that impenetrable bound- 
ary conditions at the top and bottom of the spherical 
shell lead to frequencies that are purely real (e.g., os- 
cillating waves only, with no spuriously growing modes) 
but now the eigenfunctions will be severely distorted in 
the ANS equations as compared with the LBR equations. 
In a sense, the eigenfunctions try to diverge to infinity 
but are constrained by the boundary conditions. Ana- 
lytic eigenfunctions can be found if we consider a sim- 
plified atmosphere with constant gravity g = —gr and 
constant Brunt- Vaisala frequency N . In an isothermal 
atmosphere with temperature To this can be achieved by 
setting the entropy gradient to 

V.So = ^ = f (105) 
or T 

The background entropy So is found by integration, with 
the arbitrary constant set by a reference value within the 
atmosphere (here at the base of the domain) . The back- 
ground pressure and density are determined by hydro- 
static balance and their values at the reference layer. We 
first derive the analytic solutions and then compare these 
solutions with fully nonlinear calculations using two ver- 
sions of the anelastic spherical harmonic (ASH) code. 

5.1. Modes in stratified isothermal spherical shells 

We begin by obtaining analytic solutions for the low- 
Mach number ANS, LBR, and RG equations. Full details 
are given in Appendix [Bj Motivated by the properties of 
the solar radiative zone, we solve for the eigenvalues of 
the low-Mach number anelastic equations within a spher- 
ical shell stretching from a = 0.50i?© to b = 0.70i?© with 
r the solar radius. This shell has geometric extent 

X = t = 0.717 (106) 



Table 4 

Atmosphere parameters 



n p 


AS/cp 


H 


H 


T 


N 


T"BV 






Mm 


Rq 


10 6 K 


io- 3 ^- 1 


s 


0.25 


0.1 


548 


0.788 


39.3408 


0.84 


1185 


1.0 


0.1 


137 


0.197 


9.83520 


1.69 


592.6 


2.5 


1 


54.8 


0.0788 


3.93408 


2.66 


374.6 


5.0 


2 


27.4 


0.0394 


1.96704 


3.77 


264.9 


7.5 


3 


18.3 


0.0263 


1.31136 


4.62 


216.3 


10.0 


1 


13.7 


0.0197 


0.98352 


5.34 


187.4 


12.5 


5 


11.0 


0.0158 


0.78682 


5.96 


168.5 



Note. — Quoted are the number of density scale heights 
in the domain n p , the non-dimensional entropy drop across the 
shell AS/cp, the physical size of the density scale height H 
in megameters and relative to the solar radius, the isothermal 
temperature To, the constant Brunt- Vaisala frequency N and 
the corresponding timescale tbv = 1/N. In all models, r^ot = 
a = 210A/m ss 0.30.R© and r top = b = 485Mm « 0.70i? Q , 
with x = a /b = 0.433 and with the solar radius tq = 695Afm. 



Additional simulations conducted at AS/cp 
10 — 4 are not shown here. 



10" 2 , 10~ 3 and 



Table 5 

Isothermal atmosphere solutions 



/>'4 



n p = 2.5 

ANS 9.72953 
LBR 10.1839 



ANS 
LBR 



10.4846 
12.0834 



7.5 



LBR 
ANS 
FC-1 
FC-00 



14.6998 
11.6324 
2.16693 
9.49926 



n 



P — 
ANS 
LBR 



12.5 



14.692 
20.888 



19.0635 
19.3017 



19.4637 
20.3817 



22.0671 
20.1134 
15.4867 
18.6646 



22.0618 
26.7604 



28.4829 
28.6431 



28.7527 
29.3825 



30.5764 
29.1969 
26.2203 
27.9660 



30.5726 
34.1169 



37.9246 
38.045 



38.1277 
38.605 



39.5211 
38.4639 
36.2565 
37.3421 



39.5182 
42.3167 



47.3752 
47.4717 



47.5379 
47.9216 



48.6627 
47.8080 
46.0510 
46.7676 



48.6602 
50.9577 



Note. — Radial wavenumbers for gravity waves in 
selected bounded isothermal atmospheres listed in Ta- 
ble [4] Quoted are the five lowest wavenumbers fci-fcs in 
each equation set. In the full Euler equations, the radial 
wavenumber depends on spherical harmonic l\ as such, we 
quote wavenumbers at low and high values of I (FC-1 at 
i = 1 and FC-00 at i = 50 respectively). The RG equa- 
tions have the same eigenfunctions and radial wavenumbers 
as the LBR equations and are not separately quoted. 

and we consider several different values for the scale 
height H and number of density scale heights n p . The at- 
mospheric parameters are reported in Table HI The first 
five such wavenumbers for the ANS and LBR equations 
are presented in Table [5] for several of these atmospheres. 

We begin by discussing eigenfunctions in the n p = 7.5 
atmosphere, as this atmosphere will form the primary 
comparison case for the 3-D numerical simulations in Sec- 
tion 15.21 In Figure [5] we show both the fundamental k\ 
mode and a higher-order £5 mode. In addition to the 
various low-Mach number anelastic eigenfunctions, here 
we also overplot eigenfunctions for the fully compressible 
Euler (FC) equations; and these require numeric solu- 
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Figure 2. Eigenfunctions for the n p = 7.5 atmosphere, (a) Eigen- 
functions for the fundamental fci mode and (b) the higher-order k§ 
mode, (c) Dispersion relationship u)/N for the first, third and fifth 
radial modes (fcj, kg, fcs), with lower-k having higher u). In each 
plot, the ANS equations are shown in blue (dash-dotted) while the 
LBR equations are shown in red (solid). The full compressible 
results arc shown in black. For the eigenfunctions, the solid line 
corresponds to I = 1 (FC-1) and the thick dashed line correspond- 
ing to i = 50 (FC-oo). 



tions. In the full FC equations, the radial eigenf unction 
depends on spherical harmonic £, whereas in the anelas- 
tic equation sets this coupling disappears. This effect 
is most pronounced in the FC eigenfunctions at low-£, 
with the eigenfunctions largely becoming constant with 
£ when I ^> k. As such, we plot two FC eigenfunctions, 
one at 1 = 1 (FC-1) and one at t = 50 (FC-oo). 

As is clearly evident in Figure [2j>, the discrepancies in 
the ANS eigenfunctions do not diminish at high radial 
wavenumbers. This continues to hold true for higher 
wavenumbers than we show here. This is not surpris- 
ing, as these discrepancies arise from the energy non- 
conservation in the ANS equations, rather than from as- 
sumptions about the relative size of the gravity wave- 
lengths and scale heights in the atmosphere. In con- 
trast, at high-fc, the other equation sets all converge. 
The dispersion relationship for odd modes k\, k% and 
k§ are shown as a function of spherical harmonic £ in 
Figure [2jc). For the k\ mode and at low-£, the low-Mach 



number anelastic equations generally produce higher fre- 
quencies than the full compressible Eulcr equations. At 
higher-£ all of these frequencies converge to the Brunt- 
Vaisala, frequency N, and the frequencies in the ANS and 
LBR equations generally cross the frequencies of the Eu- 
ler equations at some moderate I. The frequencies con- 
verge much sooner at high radial order (e.g., k§). The 
RG equations are not shown in Figure [2j their frequen- 
cies are consistently a factor of y^y ps 1.3 larger than the 
LBR equations. 

The eigenfunctions of the fundamental mode k\ are 
shown in Figure [3] for several isothermal atmospheres 
from Table |U With the normalization that we have cho- 
sen (Appendix [B| , the ANS eigenfunctions are generi- 
cally larger in amplitude than the other systems of equa- 
tions. This difference is most pronounced near the top of 
the domain, and the discrepancies grow as the amount 
of stratification grows. 

5.2. Numerical models with the ASH code 

We turn now to fully nonlinear 3-D simulations of grav- 
ity wave propagation using the ASH code. We study 
gravity waves in ASH using both the standard ANS equa- 
tions as well as an implementation of the LBR equations. 
In the ASH-ANS equations, the momentum and energy 
equations are 



Po 
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where the viscous stress tensor is 
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-2p is 



ij - g(V-ti)$y 



(108) 



(109) 



with eij the strain rate tensor and 5ij the Kroncckcr 
delta. These anelastic equations assume a linearized 
equation of state (eq Hip and the anelastic constraint 
(eq[5]) but are otherwise fully nonlinear. The ASH-LBR 
equations arc identical except for the momentum equa- 
tion, where 



Po 



du 
~dt 



V« = -p V — - po— g -V T>. 

\PoJ c p 

(110) 

All other properties of the simulations are identical. 

We take the geometry and atmosphere used previously 
in this section for the background reference state entropy 
So, pressure Po, temperature To, and density po. These 
quantities vary in radius but do not evolve in time. Here 
we first focus on simulations conducted in an isothermal 
atmosphere with n p = 7.5 and with other parameters 
given in Table [4] In comparison, over the same range 



of radii in the Sun n 



while n 



P,0 



6.6 over 



the whole solar radiative zone. As such, the results pre- 
sented here are likely an over-estimate for comparable ef- 
fects in the solar interior, but the larger number of scale 
heights more clearly emphasizes the differences between 
the ASH-ANS and ASH-LBR equations. Wc will return 
to solar conditions at the end of this section. 




Figure 3. Eigcnfunctions and dispersion relationships for selected isothermal atmospheres. In (a, b) n p = 2.5, in (c,d) n p = 5, and in 
(e, /) rip = 12.5. Eigcnfunctions for the fundamental k\ mode are shown in (a, c, d) and dispersion relationships for the fei , k-j, and fcs radial 
modes are shown in (6, e, /). Labels for lines in all plots are given in (e). 



In the pseudo-energy conserving ANS equations, the 
scaled pseudo-density (cq (SHJ) is weighted by the back- 
ground entropy So- As such, the non-dimensional en- 
tropy drop across the domain 

AS/cp = -L(s*o(r t op) - Soil-bat)), (HI) 

or the number of pseudo-density scale heights rip with 

rip = n p + AS/cp, (112) 

are both likely better measurements of how much pseudo- 
energies differ from energies in the ANS equations than 
the number of density scale heights n p in the atmosphere. 
For a non-isothermal atmosphere, we can use the equa- 
tion of state ((6|) to obtain 

rip = — — -n p n Ti (113) 

7 7 

where nx = In (To.bot/7o,top) is the number of temper- 
ature scale heights. In an isothermal atmosphere with 
7 = 5/3, this reduces to rip = (7/5)n p and the number 
of ANS pseudo-density scale heights always exceeds the 



number of density scale heights. Our isothermal atmo- 
sphere with rip = 7.5 has AS/c p = 3 and rip = 10.5. 

The numerical simulations were conducted in a non- 
rotating system with viscosity v = 1 x 10 10 cm 2 /s and 
entropy diffusivity k — 4 x 10 10 cm 2 /s and with cp — 
3.4x 10 8 ergsg -1 K7 1 . In contrast to ASH simulations of 
stella r convection (e.g., iBrown et al.ll2008t iMiesch et al.l 
120081) , in these isothermal atmosphere simulations we ne- 
glect radiative diffusion of temperature in the entropy 
equation and diffusion of the background entropy gra- 
dient VS*o, which in these simulations is set by equa- 
tion (|105[) . As such, there is no energy flux through the 
simulation. The velocity boundary conditions at the top 
and bottom of the domain are stress-free and impene- 
trable, and the thermal boundaries maintain a constant 
entropy gradient. All simulations are conducted with 
a resolution of 257 x 256 x 512 (N r x N e x A^), with 
all functions expanded in Chebyshev polynomials radi- 
ally and spherical harmonics horizontally; this leads to 
a dealiased spectral resolution of i max = 170, which re- 
solves the wave motions studied here. 
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Figure 4. Fluctuating velocities in numerical simulations for the 
n p = 7.5 isothermal atmosphere. Shown are the rms radial ve- 
locities for ASH simulations calculated with ANS (blue, solid) and 
LBR (red, solid) treatments of the momentum equation. The thick 
dashed lines give the analytic k\ cigenmode for each equation set, 
normalized by the peak velocity realized in the ASH simulations. 



Timcstcpping errors can have important impacts on 
the properties of wave motions. The ASH code uses 
a second-order Adams-Bashforth/Crank-Nicolson tech- 
nique for time evolution, which treats diffusive processes 
implicitly and advective processes explicitly. In our stud- 
ies here, we found that it is crucial that the advective in- 
teractions between the wave motions and the background 
reference state stratification be handled implicitly (on 
the Crank-Nicolson side). In the entropy equation (|108[) . 
this term is 

u-VSq. 

If these interactions are handled explicitly (via the 
Adams-Bashforth portion) then the solutions are sensi- 
tive to the size of the timestep; with sufficiently small 
timesteps a solution can be time-evolved correctly, but 
these timesteps must be nearly an order of magnitude 
smaller than arc otherwise possible. Larger timesteps 
lead to explicit timestep errors that grow quickly in the 
solution. Handling these interactions implicitly, as we do 
here, leads to much more stable behavior. To simplify 
matters, in these studies we fix the timestep at slightly 
less than one third of the Brunt- Vaisala timescale tbv 
(e.g., 70s in the n p = 7.5 atmosphere). 

At the start of each simulation, random entropy per- 
turbations are introduced in a band of spherical harmonic 
£ ranging from £ = 1-30, at all spherical harmonic m val- 
ues. The radial perturbation has two bumps in radius, 
defined by 



f(r) = 1 - 3a; 2 + 3a; 4 - x 6 + 2.5x - 2.5a; 3 



(114) 



with scaled radius x given by 

x = (2r - r top - r hot )/(r t0 p - r hot ), x E [-1, 1]. (115) 

This radial perturbation does not exactly match the ra- 
dial eigenfunction of the gravity waves but rather drives 
a broad band of such waves with the largest power in the 
lowest fci and ki modes. The initial perturbations lead 
to flows of roughly 1ms -1 in amplitude, with Reynolds 
numbers Re = uL/v of about 100. The viscous Q of 



Q 



uL 2 



.67 x 10 7 



v 



(116) 



where we have taken u> ~ 
1 fundamental mode has u> 
rbot , the depth of our shell 



N (the low-frequency £ 
« 0.12N) and L = r top 
A thermal Q would be four 



times smaller. The Q calculated in equation (|116p is most 
applicable to our longest wavelength modes; our shortest 
wavelength modes with ^ = 30 would have a Q of about 

Qao = * 9.32 x 10 4 , (117) 

which is still quite large. Thus, we expect that the grav- 
ity waves should only be very weakly damped by diffu- 



5.3. Eigenf unctions and violation of energy 
conservation 

We expect that there will be two clear effects from the 
choice of ANS or LBR equations. The first such effect 
is that radial eigenfunctions of the two systems should 
differ strongly, as discussed in Section 15.11 The radial 
eigenfunctions for the ASH- ANS and ASH-LBR simula- 
tions are shown in Figure [4] Plotted on the same scale 
and against radius are fluctuating rms radial velocities 
Vr' at a time late in the simulations (30 days after ini- 
tiation, or roughly 12,000 tbv)- These rms velocities 
are further time- averaged over roughly 2.5 days or about 
1000 tbv- 

Overplotted on each simulation is the appropriate ra- 
dial eigenfunction corresponding to the gravest k\ mode 
for the ANS or LBR equations. Here the eigenfunctions 
are scaled by the peak rms velocity. In both simulations, 
the rms velocities from the fully nonlinear 3-D numeri- 
cal simulations agree very well with the analytic eigen- 
functions. In the ASH-ANS set of equations, the radial 
velocities peak more strongly in the upper portion of the 
domain, reaching amplitudes 2-4 times larger than the 
ANS-LBR equations. The fluctuating velocities differ 
strongly, as expected. 

The second effect is that, as discussed in Section [4] the 
LBR equations should conserve energy while the ANS 
equations conserve a pseudo-energy. In the simulations 
we define volume-averaged total energy E, kinetic energy 
K and potential energy U densities 



K =Vj 2 Pou2dV > 

TT 1 f 1 f d S 
V J 2 \orcp 



?1) 

cp 



(118) 
(119) 

(120) 



with fluctuating velocity u and fluctuating entropy Si, 
and where the integral is over the full simulation vol- 
ume V (e.g., eg n 155)1 . Likewise we define pseudo-energy 
densities (e.g., egn lW)) 



PE = PK 

PK= — 
V 



l PU > 

e ~ So/cp Po u 2 dV, 



PU = 
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V 



(121) 
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Figure 5. Temporal evolution of energies and pseudo-energies, shown over identi cal interv als for ASH-ANS and ASH-LBR. (a) Energy and 
(b) pseudo-energy in ASH-ANS simulation, with definitions as given in equations {Tl8jfl23}. Pseudo-energy densities are here multiplied by 
10 s . The total energy E and pseudo energy PE is divided by 2 to highlight the fluctuations between K and U or PK and PU respectively. 
The ASH-ANS equations clearly do not con serve energy E but clearly do conserve pseudo-energy PE. Over the interval shown, AE 0.101 
while APE < 10 -5 , as defined in equation J 1241 ) . (c) Energy and (d) pseudo-energy for ASH-LBR simulation. Energy is clearly conserved 
in this system, while pseudo-energy is not, with AE < 10~ 6 while APE w 0.03. The temporal interval shown in all plots spans about 
lOOrgy and begins about 24,000Tgy after the start of the simulations. 



In an isothermal atmosphere with constant Brunt- 
Vaisala frequency N, So is a function of radius and the 
stratification term cannot be factored out of the inte- 
gral. If thermal and viscous diffusion can be neglected, 
E should be conserved in the LBR equations while PE 
will vary in time. Likewise, the ANS equations should 
conserve total pseudo-energy PE but should fail to con- 
serve total energy E. Indeed, this is what we find. 

Shown in Figure El are temporal traces of energy and 
pseudo-energy in the ASH-ANS simulation and the simu- 
lation using the ASH-LBR equations. Here a short inter- 
val, spanning about IOOtbvi is shown from a much longer 
simulation. The wave periods are generally longer than 
tbvi owing to their long horizontal wavelengths. The 
ASH-ANS simulation shows large variations in kinetic 
and potential energies K and U and does not conserve 
total energy E (Fig. Eta)). This simulation does how- 
ever clearly conserve total pseudo-energy PE (Fig. EK&)). 
Over much longer intervals of time, both the total energy 
and total pseudo-energy decay dissipatively. In contrast, 
the ASH-LBR simulation correctly conserves energy E 
(Fig. Etc)) while the pseudo-energy PE fluctuates in time 

(Fig. Eld)) 

We define the relative energy variation AE and 



pseudo-energy variation APE as 

AE = (5E)/(E), and APE = (5PE)/(PE), (124) 

with S signifying the standard deviation in time and 
brackets denoting a time average over the same period. 
Subtracting off the slow diffusive decay, we find that over 
a ten-day interval (4000 r B v), 5.6% in the ASH- 

ANS simulation (during the interval shown in Figure El 
AE w 10.1%). Over the same interval, APE < 0.001%. 
In contrast, the ASH-LBR simulation has AE < 0.0001% 
and APE w 3% over the same ten-day interval of time. 
We have conducted similar simulations with diffusivitics 
v and k ten times larger (e.g., Re w 10) and find a simi- 
lar level of variability, and thus conclude that our results 
are not strongly dependent on the level of diffusivity em- 
ployed. For the linear waves considered here, we find that 
AE and APE are independent of the initial perturbation 
amplitude. 

In these many-wave simulations, the introduced waves 
span varying portions of the frequency dispersion re- 
lationship, including regions where u> depends almost 
linearly on I and regions where it does not (e.g., Fig- 
ure Etc)). Thus we might expect a collection of these 
very linear waves to behave as incoherent oscillators and 
that the relative energy variations for many waves might 
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Figure 6. Relative energy variation AE in the single-wave ASH- 
ANS solution for the atmosphere with n p = 7.5 and AS/cp = 
3. This should be compared with Figure Ufa). Generally, the 
variations in these single wave solutions are about a factor of 5 
larger than the many-wave solutions. 

be smaller than those of any individual horizontal wave. 

This is confirmed by simulations where only a single 
spherical harmonic perturbation is initially introduced, 
as shown by traces of E, K and U in Figure [S] for the 
rip = 7.5 and AS/cp = 3 atmosphere. In this ASH- 
ANS simulation, only £ = 30 waves (at all \m\ < £) are 
initially excited, with the same radial perturbation as 
the many- wave simulations (eq 1 1 14[) . Hereon, we will 
refer to these as single-wave solutions. Comparing Fig- 
ure [5] with the corresponding many- wave solution in Fig- 
ure [5ja) it is clear that the relative energy variation is 
significantly greater. Here, AE pa 25% (APE remains 
negligible). We have studied single- wave solutions with 
different horizontal wavelengths, sampling in the range 
from I = 1 to £ = 100 and find that this level of en- 
ergy variation is reasonably representative for individual 
waves of any horizontal wavelength in this range. This 
confirms our understanding that the phenomenon of en- 
ergy non-conservation is due to the level of stratification 
and docs not depend strongly on the particulars of any 
single mode (e.g., horizontal or vertical wavelength). 

We now turn to considering isothermal atmospheres 
with differing levels of stratification, ranging from n p = 
0.1-12.5 and AS/c P = 0.1-5 (see Table HJ. The configu- 
rations of the simulations are the same as previously dis- 
cussed, though at large stratification (n p > 10) a higher 
resolution was used, with N r = 1025 and a dealiased 
horizontal resolution of £ max = 340. 

The time-averaged relative energy variations in these 
ASH-ANS simulations are shown in Figure [3 which 
displays both many- wave solutions (£ = 1-30, tri- 
angles) and single-wave solutions (£ = 30, squares). 
The many- wave solutions span from AS/cp = 0.4-5, 
while the single- wave solutions span a wider range from 
AS/cp = 0.0001-5. As the stratification increases, en- 
ergy non-conservation becomes increasingly significant in 
the ASH-ANS equations, with AE approaching 10% in 
the many-wave simulation with n p = 12.5 and AS/cp = 
5 and AE sa 43% in the corresponding single-wave solu- 
tion. Generally, we find that the relative energy varia- 
tions are about 5 times higher in the single- wave solutions 
than the corresponding many-wave solutions, indepen- 



dent of stratification. In the corresponding ASH-LBR 
simulations (not shown), energy is always well conserved 
with AE < 10- 6 . 

Unsurprisingly, the relative energy variation is smaller 
in less stratified atmospheres. At very low levels of 
stratification (AS/cp — > 0) the ratio of pseudo-density 
and density is almost constant throughout the domain 
feq [89|) , and we should expect the pseudo-energy conserv- 
ing ASH-ANS equations to also conserve energy fairly 
well. Indeed, this is what we find. As shown for single- 
wave solutions in Figure EJa), at low values of AS/cp, 
the energy variation AE is also small (AE pa 5 x 10~ 5 at 
AS/cp = 0.001). With increasing stratification, AE in 
the single- wave ASH-ANS solutions scales almost linearly 
with AS/cp up through AS/c P = 0.1. At AS/c P = 0.4 
AE pa 1% in the single-wave solution and AE sa 0.2% 
in the many-wave solution. There is a change in the 
scaling for both single-wave and many-wave solutions at 
AS/cp pa 1, apparent in both Figures [7Ja, 6); we do not 
understand the origin of this behavior. In all cases shown 
here, APE < AE, and generally APE ~ lO^-lO" 6 . 
This floor on APE likely reflects aspects of our data 
analysis technique and we feel that our current approach 
is insufficient to reliably measure energy variations in 
cases where AS/cp < 10 -4 . 

5.4. Nonlinear interactions 

At much lower levels of diffusivity, or at larger ini- 
tial amplitudes, the gravity waves may begin to interact 
nonlincarly. This is also likely to occur when the gravity 
waves are driven by overshooti ng convection from below 
(e.g., IMihalas fe Toomrd 119811 ). To confirm the linear 
nature of the waves we have studied here, we define a 
Froudc number Fr as 



or the ratio of local vorticity to the Brunt- Vaisala fre- 
quency N. Thi s corresponds to th e vorti city criteria for 
nonlinearity in Mihala s fc Toomrel ([19811 ) . We find here 
that Fr attains a peak value of about 5 x 10~ 5 in the 
n p = 7.5 ASH-ANS simulation and of about 1 x 10~ 5 
in the corresponding ASH-LBR simulation. Thus the 
waves studied here are quite linear. For linear waves, 
the Froude number gives the characteristic amplitude of 
all fluctuations. Owing to this, despite the large strati- 
fications studied here (AS/cp ~ 1), the thermodynamic 
fluctuations remain quite small (S±/cp ~ 10 -5 ). 

When nonlinear interactions become important, we 
might expect that the non-conservation of energy may 
cause ASH-ANS simulations to diverge even further from 
simulations which do conserve energy (e.g., the ASH- 
LBR equations). As discussed in Section |4~21 the conser- 
vation of pseudo-energy also vanishes when nonlinear- 
ity is important. Energy is conserved in the nonlinear 
LBR equations, but neither the pseudo-energy nor the 
energy is conserved in the nonlinear ANS equations. If 
the pseudo-energy is also not conserved, it may be pos- 
sible to inject pseudo-energy into otherwise closed sys- 
tems; alternatively, the energy and pseudo-energy may 
leak away without coupling to the reservoir of internal 
energy. Either case leads to physical inconsistencies. 

Lastly, the transport by nonlinear processes in the 
ANS equations is likely to be very different from that 
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Figure 7. Violation of energy conservation in ASH-ANS simu- 
lations in various atm osphe res, (a) Relative energy variation AE 
as given in equation II124H for isothermal atmospheres with dif- 
ferent non-dimensional energy drops AS/cp in log-log plot, (fe) 
Same, in linear plot, emphasizing the behavior at large AS/cp. 
Shown in both are solutions with a single horizontal wave (I = 30, 
blue squares) and solutions with many horizontal waves (I = 1-30, 
black triangles) . Also shown arc solutions in a solar radiative zone 
atmosphere stretching from 0.5-0.7-Rq, with a single-wave solu- 
tion (blue square with asterisk) and a many-wave solution (black 
diamond with asterisk). All solutions arc time-averaged over an 
interval of 2000ra^, generally beginning about 20Qrpy after the 
start of the simulation. 

in equations that do conserve energy, as the eigenfunc- 
tions of gravity waves in the ANS equations are signifi- 
cantly higher in amplitude in the upper domain of stably- 
stratified atmospheres (Figures [2HI]) ■ This will also have 
important implications for mode coupling, for the steep- 
ening and breaking of gravity waves, and for all other 
problems where the shape of the eigenfunction itself is 
important. 

5.5. Solar atmospheres 

The results presented here so far have been for the 
special case of an isothermal atmosphere. The solar ra- 
diative zone is stably stratified, but has fewer density 
scale heights than have been considered in most of this 
section. Across the entire solar radiative zone 

M.7Rq -i Qo 

(AS/c p ) & = / __p«i.3, (126) 

and rip,© ps 7.9, while AS/c p ps 0.42 and np t Q ps 2.3 over 
the shell geometry that we consider here (0.5-0.7i? Q ). 

To constrain our results, we have repeated these grav- 
ity wave rundown ASH-ANS and ASH-LBR simulations 



in the solar radiative interior. We take our mod el at- 
mosphere from the CESAM code (|Brun et al.ll2002fl . We 
keep the same values of v and k and continue to ne- 
glect radiative diffusion acting on the fluctuating flows. 
A large scale radiative diffusion based on the Rossland 
mean-opacity is included that acts on VTo and there is a 
flux equal to the solar flux throughout the domain. We 
keep the same choice of r top and rbot, thus n p ps 1.8 and 
AS/cp ps 0.417. 

These solar simulations are shown in Figure [7ja) as 
asterisks. In this model solar atmosphere, we find that 
in many-wave solutions (I = 1-30) with the ASH-ANS 
equations AE ps 0.9% and APE is tiny while in the 
ASH-LBR equations APE ps 0.9% and AE is tiny. 
In single- wave solutions (I — 30), the relative energy 
variations in the ASH-ANS equations are much larger 
(AE ps 4.5%). Surprisingly, the solar atmosphere simula- 
tions show larger relative energy variations than similarly 
stratified isothermal atmosphere simulations, with AE 
being roughly five times larger in this solar atmosphere 
than in the corresponding AS/cp = 0.4 isothermal at- 
mosphere. If we plotted these against rip, the solar simu- 
lations would lie midway between the AS/cp = 0.1 and 
0.4 atmospheres and would still be clearly discrepant. 
We expect that the effects of energy non-conservation 
will become significantly larger as more entropy scale 
heights arc included in the domain. This may be dif- 
ficult to diagnose in simulations that include a realistic 
solar stratification as the radially varying Brunt- Vaisala 
frequency creates acoustic cavities that may trap high 
frequency gravity waves, but we expect that the low fre- 
quency waves which travel the entire radiative zone and 
experience the full stratification will be affected. 

6. RECOMMENDATIONS FOR IMPROVING 
ANELASTIC TREATMENTS OF STELLAR 
INTERIORS 

The results of Sections [5H5] provide a clear path to 
improving the treatment of dynamics within stably- 
stratified atmospheres in anelastic systems of equations. 
As clearly shown in Figure El the ANS equations do 
not conserve energy and instead conserve a stratification 
weighted pseudo-energy. These equations thereby obtain 
incorrect frequencies and radial eigenfunctions for grav- 
ity waves in both infinite and bounded isothermal atmo- 
spheres. The variation in the eigenfunctions is substan- 
tially larger than the level of energy non-conservation. 
These results hold in general for all subadiabatically- 
stratified atmospheres. 

In contrast, the anelastic LBR equations do conserve 
energy and appear to need no additional modification to 
capture dynamics in subadiabatically-stratificd regions. 
This is fairly surprising, as those equations are generally 
derived in nearly adiabatic atmospheres and the isother- 
mal atmospheres we have considered here take them far 
from their realm of validity. At low vertical and horizon- 
tal wave numbers, eigenfunctions in the LBR equations 
differ from the full compressible equations, and results 
from gravity waves in this regime should be treated with 
caution. These differences shrink as either wavenumber 
increases and the LBR equations may do a reasonable job 
of capturing the dynamics of those shorter wavelength 
gravity waves (e.g., Figure [2]). Though we did not ex- 
plore their dynamics in direct numerical simulation, we 
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have demonstrated that the RG equations do not con- 
serve energy in general atmospheres that have tempera- 
ture gradients, and in those atmospheres will also likely 
obtain incorrect radial eigenfunctions. 

To correctly capture the dynamics in sub-adiabatically 
stratified regions, it is vitally important that subsonic 
treatments of the fluid equations conserve energy. Sys- 
tems of equations that conserve energy are physically 
self-consistent, even if simulations done with them have 
transport coefficients (e.g., v and k) that are several 
orders of magnitude larger than the molecular values 
in astrophysical systems. Systems that do not con- 
serve energy are not physically consistent, and though 
the variations of energy may be small for some prob- 
lems, these variations point to deeper underlying prob- 
lems with those systems of equations. In particular, the 
eigenfunctions of the waves are significantly different in 
the non-conscrvativc systems (e.g., ANS) from the energy 
conserving systems (e.g., FC and LBR), and this is very 
important for nonlinear transport, mode coupling and 
wave steepening and breaking. The clear path forward 
is to ensure that simulations employ anelastic systems of 
equations that conserve energy; fortunately this can be 
done with simple modifications to existing codes. 

The route to energy conservation is to modify the mo- 
mentum equation of the non-conservative systems. Fun- 
damentally the conservation of energy is more physical 
than the conservation of momentum: there are many 
physical systems conserve energy instead of momentum, 
especially those where a very fast restoring force acts to 
constrain the behavior of the system. Examples of this 
include inelastic scattering off of rigid boundaries, where 
momentum changes sign but energy is conserved, and 
roller coasters on rigid tracks, where the track changes 
the momentum of the careening roller coaster but not 
its total energy. In anelastic systems, the fast sound 
waves provide the rapid restoring force and the diver- 
gence constraint embodied in equation © acts analo- 
gously to the rigid tracks of the roller coaster, applying 
a continuous forcing to the system. Fundamentally, this 
forcing is energy-conserving in the LBR equations but 
violates energy conservation in the ANS equations. 

In simulation codes, there are two equivalent paths to 
convert the ANS equations into an energy conserving 
form identical to the LBR equations. The first path is 
by rewriting the equations to exactly match the anelastic 
LBR equations. This is done by solving for the reduced 
pressure w instead of the fluctuating pressure Pi, and 
by converting the buoyancy term to a "codensity" where 
entropy fluctuations S\ contribute to buoyancy but pres- 
sure fluctuations do not. Doing so causes the momentum 
equation to take the following form 



du _ 
— + u ■ Vu 

at 



-Vzu g 



(127) 



which is identical to the LBR momentum equation (| 1 10[) . 

The second path to energy conservation is consider- 
ably simpler and relies on introducing a correction term 
into the momentum equation. Generally, energy non- 
conservation occurs in anelastic systems of equations 
when the reduced pressure w interacts with the back- 
ground stratification. The problematic term in the ANS 
momentum equation (p~8j) is the zu'V(So/c p ) term. We 



modify the momentum equation to read 
^+wV« = -— VPi + f + P B vz ) g-—g-VT>, 



dt 



Po 



7P> 



where the correction term is 

P 



Fbvz = —V (So/cp) . 

gpo 



(128) 
(129) 



and where equation (|128|) reduces to the LBR momen- 
tum equation (|110[) . We remind the reader that we have 
taken g = —gf, and this sign is incorporated into equa- 
tion (|129|) . In a code like the ASH code, where the in- 
termediate variable p\ is carried around, this amounts to 
changing the equation of state for density fluctuations to 



Pi 
Po 



Pi 



Pi 



gpo 



Si 



Si 



= Vlnpo - 

cp J c P gp c P 

(130) 

Equations (|127|) and (|128p arc mathematically equiva- 
lent, but implementing this second path in a production 
code like ASH requires only a few lines of code and is con- 
siderably simpler than re- writing the equations in terms 
of za. We have implemented both approaches in the ASH 
code, and the two paths to energy conservation give iden- 
tical results in these test simulations. 

The RG equations appear to not conserve energy in 
any atmosphere that has a temperature gradient. Instead 
they conserve a pseudo-energy weighted inversely by the 
background temperature To in both stably-stratified ra- 
diative zones and in nearly adiabatically-stratified con- 
vection zones. This is an important distinction, as the 
differences between pseudo-energy conservation and en- 
ergy conservation can appear in the RG equations even 
when AS/ cp is small. Instead, what matters is the num- 
ber of temperature scale heights across the domain. This 
can be seen from the RG pseudo-density (|103[) and the 
associated pseudo-density scale height 



Si 



n p - n T , 



(131) 



where hit is the number of temperature scale heights, 
(compare with cqn I113| ). In the RG equations, tit has 
a similar but opposite role as AS/cp in the ANS equa- 
tions: increased fir leads to less pseudo-density strat- 
ification. In the solar interior, To is about 15 x 10 6 K 
near the core, about 2 x 10 6 K at the base of the con- 
vection zone (0.7i?o) and roughly 4 x 10 5 K in the upper 
convection zone (0.93Rq) . This leads to ut ~ 2 and 
rip « 4.7 across the solar radiative zone (0. 001-0. 7Rq), 
and to tlt ~ 1.6 and rip w 0.9 across the deep convection 
zone (0.7-0.93i? Q ). The energy conserving properties of 
the RG equations could be studied in either a solar inte- 
rior setting or in a polytropic atmo sphere where ther e is 
a linear temperature gradient (e.g.. lJones et al.ll2011| ). 

As with the ANS equations, it is straight-forward 
to put the RG equations into energy-conserving form. 
The term that leads to energy non-conservation is the 
zuV In To term in equation (f2"B|) . which arises from a cor- 
rection term, 



1 Pi dT . 

7^ ^-gr = zuVlnTo, 

T gpo or 



(132) 



intended t o more correctly capture sub -adiabatic strat- 
ifications (jRogers fc Glatzmaierll2005bD . If the RG mo- 
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mcntum equation were re-written as 

^+u- Vu = -V (H) +^gr, (133) 

then these equations would conserve energy, though the 
frequencies of gravity waves in these equations remain 
a factor of y 77 / higher than both the LBR frequencies 
and the low-frequency branch of the full compressible 
Euler equations. The source of this remaining disagree- 
ment remains unclear. It would be very interesting to see 
whether the non-conservation of energy has any impacts 
on the nature of convection in the RG equations, and 
on the coupling of convection to stably-stratified regions 
above and below. 

Here we have explored how gravity waves in a solar ra- 
diative interior may be affected by anelastic treatments. 
The Boussinesq equations, which we have not consid- 
ered here, are well known to conserve energy both lin- 
early and nonlinearly, and this indicates that the issue 
is not the filtering of sound waves alone. Rather, it is 
the treatment of filtered, subsonic motions in a stratified 
atmosphere (in Boussinesq treatments the background 
density is constant). Other treatments of subsonic mo- 
tions, such as the pseudo-incompressible equations and 
Reduced Sound Speed Techniques (e.g.. iRempell [20051 
120061 : iHotta et al.l 120121) may similarly not conserve en- 
ergy, and we would suggest that this be carefully tested. 
Generally, the isothermal atmospheres considered here 
or similar stably-stratified polytropic atmospheres (not 
considered here) provide simple test cases. The next pa- 
per in this series will consider variations on the pseudo- 
incompressible equations. 

The subsonic dynamics of gravity waves may play an 
important role within the radiative envelopes in more 
massive stars, such as main-sequence A-, B- and O-type 
stars, where convective overshoot drives gravity waves 
up into a rarifying envelope leading to possible nonlinear 
wave breaking. Taking a CESAM model of a 2M Q A- 
type star, we estimate that the entropy change across 
the radiative envelope is about 

MS7R, i an 

(AS/c P ) 2Me = / — ^«3.8, (134) 
Jo.i5R, c P dr 

which is substantially larger than the drop across the so- 
lar radiative zone (eq II 26[) . Over this range of radii, 
n p w 16 and rip w 20. Though we have focused 
on stellar interiors, energy conservation within anelas- 
tic systems is an important concern for dynamics in any 
stably-stratified atmosphere, including planetary atmo- 
spheres, pl anetary interiors, and astr ophysical accretion 
disks (e.g.. iBarranco fc MarcusH 2005) . Energy conserva- 
tion in anelastic equations may also play an important 
role when magnetic fields arc included in questions of the 
dynamics of magnetohydrodynamic instabilities includ - 
ing magnetic buoyancy instabilities (Berkoff ct al.|[2010l ). 

Conservation of energy remains among the most sacro- 
sanct and useful of principles in the physicist's toolbox. 
The principle of energy conservation applies not only to 
ideal systems but furthermore to fundamentally dissipa- 
tive or externally driven situations. In the latter sce- 
narios, a time-dependent statement of energy budget re- 
places the simpler notion of time-constancy of total en- 
ergy in an isolated ideal system. In all situations, we 



believe that one should not tolerate the existence of un- 
controlled spurious kinetic sources. For the particular 
problem of gravity and acoustic waves, this principle is 
more than philosophical. When examining the gravity 
and acoustic waves, energy-conserving anelastic models 
reproduce fully compressible results with much greater 
fidelity than energy-violating anelastic models; and this 
fact produces implications for our understanding of stel- 
lar interiors. Our particular work further shows that the 
existence of a conserved "pseudo-energy" does not res- 
cue the energy non-conserving models. Rather, the ex- 
istence of the pseudo-energy merely indicates why these 
problems appear to have gone unnoticed in previous sim- 
ulations. However, including any nonlinear or dissipativc 
effects leads to the impossibility of even proper pseudo- 
energy budgeting; in this case dissipation can even inject 
positive pseudo-energy. There is in short no way around 
the issue. Our best advice: properly account for energy 
whenever possible. 
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APPENDIX 

A. CONSERVATIVE FORM OF BUOYANCY 
TERM 

When analyzing energy conservation properties in the 
anelastic equations (Section|4]), the energy equation takes 
the general form 

S 

— + V • [u(K + Pq w) }+p u-g— = RHS, (Al) 
ot c p 

with kinetic energy density K = pou 2 /2 and where 
RHS is t he r ight hand side (e.g., cqns [Ml and l93| . In 
equation (|A1[) . the buoyancy work term takes the form 
poU-gSi. We put this buoyancy work term into conserva- 
tive form by using the gravitational potential, g = — V$ 
and by recognizing the relation 

p u ■ gSi = -d t (p $Si) - V • (poM^Si) - $ Pou ■ VS . 

m _ (A2) 

Using equation (|A2[) . the left hand side fcqn. IA1| ) can be 
put into conservative form 
dE 

— + V-[u(E + p m)]=RHS (A3) 
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where 



E 



Pa 



2 



w = w-— f ${r')dS {r'), 



(A4) 
(A5) 



with a an arbitrary reference radius (here the radius of 
the lower boundary) and where the entropy profile is 
monotonic. 

The energy defined in equation (|A4|) is slightly unusual 
in that the buoyancy contribution is not quadratic in 
entropy perturbation Si. This can be put in a more 
familiar quadratic form by first noting that for arbitrary 
(possibly nonlinear) motions 

Si rf$ Si „_ po d<S> DSf 
pau-g — = -p — u • Vb - 



dS c p 



2c„ dSg Dt 



d t A + V ■ (u A) -Au- Vln[ — ) , (A6) 



where the available potential energy A is given by 



and 



A = 1 po d$ s2 
2 c p dSo 1 ' 



d$ g(r) 2 



dS c p N{rY ■ 
Equations (|A6HA8[) let us rewrite equation (|A1|) as 



(A7) 
(A8) 



dE 
~dt 



V • 



u(E + pqw) 



Au ■ V In 



dS 



where the alternative total energy E is 



E= -po\u\ 



Po_d®_ q2 

2c p dS 1! 



RHS, 
(A9) 

(A10) 



and where vo is given by equation ()A5|) . 

The left-hand side of equation (jA9| can be put into 
conservative form if the condition 

/ d<S> 



Au • V In 



\dS 







(All) 



is satisfied. This happens under under two different con- 
ditions: (i) if d<f>/dSo is constant, or (ii) we only consider 
linear perturbations. If condition (i) is satisfied (e.g., 
in isothermally- or adiabatically-stratificd atmospheres) 
then equation (jAlip holds for nonlinear motions as well 
and systems of equations with RHS = will conserve a 
quadratic potential energy for nonlinear as well as linear 
motions. 

B. EIGENFUNCTIONS FOR A BOUNDED 
ATMOSPHERE 

Our analytic approach is similar to that in an infi- 
nite isothermal atmosphere, except now the wavelike per- 
turbations are expanded in spherical harmonics and the 
radial cigenfunctions must be solved for. We take the 
spherical shell geometry of Section [5] and take impenetra- 
ble boundary conditions at the upper and lower bound- 
ary 

£ r = at r = a, b (Bl) 



where a — rbot and b = r top (and see Table 3]). Analytic 
cigenfunctions can be found if we consider a simplified 
atmosphere with constant gravity g = —gf and constant 
Brunt- Vaisala frequency N, and we do so here as well as 
in the main body of the text. 

B.l. LBR eigenf unctions 

We begin with the LBR equations. In this system, 
in a spherical shell geometry, equation (|52|) for reduced 
pressure w becomes 



1(1+1) 



(B2) 



where we have used the anelastic continuity equation 
Defining 

<P(r) = t r (r)re { - r/2H) 



and 

X = -£(£+1)(1-N 2 /oj 2 ) 
the momentum equation (|60|) becomes 



d 



d 



AH 2 



(r) = A0(r). 



(B3) 
(B4) 

(B5) 



For different A the different <fi(r) arc orthogonal and this 
can be used to determine 

A = i + fc 2 , (B6) 

where the vertical wavenumber k is normalized by the 
pressure and density scale height H . With equation (|B4|) 
we obtain the dispersion relationship 



N 



(B7) 



£{£ + l) + k 2 , 
The vertical wavenumber k n can be approximated as 

,2 n 2 n 2 ( b 2 -a 2 ln(£) \ ml . 
k l = —772 1 + 0^2 I I {a \,,2 + OiH-') 



In I 



8H 2 nW+foi 



(B8) 

for a spherical shell with lower boundary at r = a and 
upper boundary at r = b. Exact solutions can be found 
numerically by solving the problem in terms of Bessel 
functions, with 



r.LBR 



(r) 



\2HJ 



\2HJ 



^ 3/2 exp(^), (B9) 



\2HJ \2HJ 
where In- and are modified Bessel functions of the 
first and second kind respectively with imaginary index 
ik n . The impenetrable boundary conditions at r = b 
requires that £(r) =0 and thus 



Kik hk ( — 



Iik (2H) Kik ^2H) _0 ' 

(BIO) 

which can be solved by Newton's method and using equa- 
tion (|B8|) as an initial guess, yielding k n . 



20 



Brown, Vasil & Zweibel 



B.2. ANS eigenf unctions 

At this point the eigenfunctions for the linear ANS 
equations can be found by a simple transformation 

f->£exp(-S /cp), (Bll) 

H -> jH, (B12) 

which leads to 

w — > tucxp (— So/cp), (B13) 

and transforms the linearized LBR wave equations (|60[) 
into the linearized ANS wave equations (ppTj) . This trans- 
formation leads to eigenfunctions of 



£r,ANsM 



a 



2-jH 
r 



In 



2 1 H 



— S /2 

r 1 exp 



(27 -1) 



2 7 i7 



(B14) 

The vertical wavenumbers are found as before by solving 



K. 



i k 



2jH 



Hk 



2 1 H 



Hk 



2 1 H 



( b 



\2jH 



(B15) 



B.3. RG eigenfunctions 

The linear RG are already in almost the same form as 
the LBR equations, except that 

X RG = + 1) (1 - 7 A 2 /^ 2 ) = \ + k 2 . (B16) 

Thus the solutions for radial wavenumber k and the ra- 
dial eigenfunctions are the same as in the LBR equations, 
but the frequencies u are a factor of y/j higher than the 
Brunt- Vaisala frequency. 

B.4. Normalization of eigenfunctions 

As defined so far, the amplitude of the eigenfunctions £ 
is a free parameter. In all sets of equations, we normalize 
the eigenfunctions by an amplitude A, with 



/ r ) 2 ex P[ — er/H]r 2 dr 
J cxp[— er / H]r 2 dr 



(B17) 



where e represents the imaginary part of the vertical 
wavenumber K and is 

'2 - (I/7) ANS equations 
1 all others. 



(B18) 



This choice of normalization gives the correct amplitude 
for motions in the different systems of equations when 
subject to the same initial conditions (entropy pertur- 
bations of fixed initial amplitude); this is how we con- 
duct the 3-D numerical simulations and thus the ana- 
lytic eigenfunctions shown in Section 15.11 show the same 
amplitude ordering as the numerical simulations of Sec- 
tion EH 
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